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Abstract. - We study a simplified, solvable model of a fully-connected metabolic network with 
constrained quenched disorder to mimic the conservation laws imposed by stoichiomctry on chem- 
ical reactions. Within a spin-glass type of approach, we show that in presence of a conserved 
metabolic pool the flux state corresponding to maximal growth is stationary independently of the 
pool size. In addition, and at odds with the case of unconstrained networks, the volume of optimal 
flux configurations remains finite, indicating that the frustration imposed by stoichiometric con- 
straints, while reducing growth capabilities, confers robustness and flexibility to the system. These 
results have a clear biological interpretation and provide a basic, fully analytical explanation to 
* ^ features recently observed in real metabolic networks. 
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Metabolic networks represent the biochemical machin- 
ery by which cells dispose of the nutrients found in their 
surrounding environment to produce the macromolcculcs 
needed for their survival, including nucleic acids, mem- 
branes, cell walls, proteins and the energy carrier ATP. 
Understanding how functionality emerges from such a 
complex system of biochemical reactions is a major is- 
sue with important implications in biocnginccring and 
pharmacology. Unfortunately, detailed kinetic models of 
genome- wide metabolic networks are computationally pro- 
hibitive and fine tuning of parameters would require an 
amount of biological information that is not available. It 
is therefore important to develop methods to predict the 
organization of metabolic fluxes from the sole knowledge 
of the stoichiomctry, which represents the full available 
information about the network topology and about the 
proportions in which different reagents are interconnected 
as a set of input-output relations. This problem is cen- 
tral in systems biology. A key question in this context is 
whether the physiological states observed in real cells are 
optimal from a metabolic viewpoint, i.e. if reaction fluxes 
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self-organize so as to maximize the production capabilities 
of the network (or of a smaller set of metabolites) [1-4]. 

It is quite intuitive that, in a generic situation where 
reagents flow in a heterogeneous input-output network, 
the ways in which fluxes can be chosen will be heavily 
constrained by the need to meet prescribed production 
objectives, while the dependence of the global properties 
on the specific realization of input and output coefficients 
will become weaker and weaker for larger networks. Ul- 
timately, it is reasonable to think that the feasibility of 
flux configurations meeting all constraints, the robustness 
of the solutions against localized flux variations (or reac- 
tion knock-outs) , and the productive efficiency of solutions 
will depend crucially on structural parameters like the ra- 
tio between the number of reactions (variables) and that 
of reagents (constraints) . 

But in order to address the question of optimality in a 
realistic biochemical network it is necessary to take into 
account the peculiar nature of stoichiometric coefficients, 
which enforce mass balance conditions at each reaction 
node in the network. Such relations in turn define a fam- 
ily of conservation laws reflecting the existence of dynam- 
ical invariants of purely topological origin that are able 
to affect the kinetics of the system as a whole [5] . To be 
explicit, let af (respectively, 6f) denote the output (re- 
spectively, input) stoichiometric coefficient of metabolite 
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/i G {1, . . . , M} in reaction i e {1, . . . , N}. A group G of 
metabolites satisfying 

£«-&?)=0 Vz (1) 
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is such that the total number of molecules in the pool, 
or the aggregate concentration, does not change in time. 
Such pools are ubiquitous in real metabolic networks: an 
important example is the total quantity of ATP and ADP 
(the 'adenylate moiety') which remain constant during 
metabolic activity, as ATP is continuously discharged to 
ADP and vice- versa [5]. 

In this Letter we study how the presence of conserved 
metabolic pools, reflecting the underlying mass-balance 
conditions of stoichiometric origin, affect the global growth 
capabilities and the volume of the solution space in a toy 
metabolic network where reagents are fully interconnected 
and stoichiometric coefficients (which play the role of a 
quenched disorder) satisfy laws such as (1). To be as gen- 
eral as possible and to highlight the ingredient of con- 
served pools, we adopt the framework of Von Neumann's 
expanding model [6], under which one aims at comput- 
ing the maximum value of p > for which the system of 
inequalities 

N 

5>«-p6f)>0 V M (2) 
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admits solution(s) in the form of non-zero flux vectors 
s = {si > 0}. From a physical point of view, the con- 
ditions (2) ensure that the total output of each metabolite 
p is at least p times the total input, so that if the maximum 
achievable p (which we shall denote by p*) exceeds (respec- 
tively, is smaller than) one, the network is in an expanding 
(respectively, contracting) state, p* = 1 signals instead a 
stationary network. (See [7] for a more thorough descrip- 
tion of Von Neumann's setup and a dynamical derivation 
of (2).) 

Following [7], we assume that every reaction produces 
and consumes (in different proportions) every metabolite. 
In such a fully-connected framework, when N — > oo and 
the input and output matrix are independent, identically 
distributed quenched random variables, it has been shown 
that p* depends only on n = N/M in such a way that 
the system undergoes a transition from a contracting to 
an expanding phase at n c = 1. Moreover, (2) admits 
a unique feasible flux configuation when p = p*. The 
advantage of using an unrealistic fully-connected setup lies 
mainly in its analytic tractability. We shall see that indeed 
the fully connected model provides an excellent proxy for 
a real metabolic network, at least as long as production 
capabilities and solution space are concerned. Graphical 
versions of the model yield essentially the same scenario 
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To keep things mathematically simple, we account for 
conserved metabolic pools by constraining each disorder 



sample (i.e., each realization of the stoichiometric coeffi- 
cients) to embed a pool formed by a finite fraction <j> of 
metabolites. In other terms, we request that 

M 

5>>f-6f) = Vz (3) 
n=i 

where is a quenched random variable that equals 1 with 
probability (f> and zero otherwise. In this way, we include 
a single conserved pool in the system. It will become clear 
that both the number of pools and their size do not affect 
the growth properties in the framework we consider. 

For a start, note that multiplying each term in (2) by 
z^ and summing over p one easily obtains 

N M N M 

5> £ z "(< - b *) + (! - p) E s « E > o ( 4 ) 
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If the z^'s are such that (3) holds, then cither p < 1 or all 
Si 's connected to a metabolite in the conserved pool must 
vanish. The latter condition however leads to the null 
solution Si = Vi in a fully connected model and must 
be discarded. Hence necessarily p < 1. One then sees 
that an expanding regime cannot occur in a system with a 
conserved metabolic pool, independently of its size. This 
suggests a radically different and more realistic scenario 
than the unconstrained case. In addition, it is easy to un- 
derstand that the inequalities (2) must become equalities 
at p — 1 for all metabolites belonging to a conserved pool. 
In other words: if a metabolite belongs to a conserved pool 
there can be no net production or consumption of it in the 
optimal state. 

To get a more thorough insight, it is necessary to com- 
pute p* exactly as a function of n. As in [7], the calcula- 
tion can be carried out in two steps. First, we compute the 
typical volume of flux configurations compatible with (2) 
and (3). This requires the calculation of an average over 
the constrained quenched disorder {af,6f}, for which we 
shall resort to the replica trick [9]. Next, following Gard- 
ner [10], we impose that the average distance between so- 
lutions vanishes (i.e. that the typical volume reduces to 
a single point). This condition intuitively marks p* , if an 
increase of p reduces the set of flux vectors satisfying (2) 
until no more solutions are found. We will see that this 
Ansatz describes the system correctly only up to a criti- 
cal value of n. Above this point, the solution is no longer 
unique. The breakdown of the Ansatz is directly linked to 
the existence of a conserved metabolic pool. 

The volume of flux vectors satisfying (2) is given by 




where the extra factor (here and in (7)) is added 
for convenience. By analogy with known systems with 
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quenched disorder, we expect the typical volume of solu- 
tions at fixed p to be given by W{p) <~ e Nv ^ , where (with 
over-bar denoting the average over the quenched disorder) 



v{p) = lim — \ogV(g) 

N— *oo iV 
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is a self-averaging quantity. Note that the disorder- 
averaging includes the constraints (3), so 
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where (• • • ) is now an average over the free, un-constraincd 
stoichiometric coefficients In turn, the quenched 

average can be computed via the replica trick: 



log V(g) = lim - log V(g) r 

r— >0 r 
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As in [7], it is convenient to write the stoichiometric 
coefficients as of = 5(1 + £f ) and 6f = 5(1 + rft ), with £f 
and zero-average Gaussian random variables. Inserting 
these in (2) one immediately sees that, to leading order in 
N, the optimal growth rate depends only on the average 
input and output coefficients: p* = a/b. In turn, the 
corrections to the leading order can be captured by re- 
writing p as 

5 /„ g 

P= r 1 + -7= 



(9) 



(3), however, requires that the average input and output 
coefficients are the same (this can be easily seen by direct 
substitution). Hence we shall set p <~ e g /^ and shift the 
focus to g: g > (resp. g < 0) now signals an expanding 
(resp. contracting) phase. The calculation in the limit 
N — ► oo can be carried out along the lines of [7], except 
that the new constraints (3) and the fact that the average 
flux is now a free variable (in [7] it was conveniently fixed 
to 1 because of the invariance of (2) under re-scalings — > 
Xsi Vi) lead to the introduction of extra order parameters. 
The key one is however still the overlap 
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between different solutions s a = {s ia } and sp = {s^} (at 
fixed g). In the replica-symmetric approximation (which 
is putatively exact in the present case due to the convexity 
of the space of solutions) , 



q a p = q + x$a 
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and one is lead to consider the following saddle-point prob- 
lem: 



v (g) = max 

m,q,x,u 



Hi (m,q,x, u) + 

- max H 2 (m,m,q,x,P,T,u,u) (12) 

p,T,m.u 



where 



H x = 



log 



dc 



2k X 



(13) 



and 
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In (13) and (14), (■■ denotes an average over the unit 
Gaussian random variable £, while k = (?yf — ) 5 



is a 



parameter related to the quenched disorder distribution. 

Physically, the order parameter x measures the dis- 
tance between solutions (it is easy to see that indeed 
(1/iV) J2i( s ia ~ Si/3) 2 = 2x(l - 5 a p)). It is reasonable to 
expect that as g increases, the typical volume of solutions 
shrinks as it gets more and more difficult to satisfy (2). 
Assuming that a single flux state satisfies all constraints 
when g — g* is then equivalent to studying the solution of 
(12) in the limit x ~~ * 0. This can be done by introducing 
a proper re-scaling of the order parameters in terms of % ; 
similarly to [7], so as to allow the integrals in (13) and 
(14) to be evaluated by steepest descent. In the present 
case, it is convenient to set 
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With these definitions, one obtains the following set of 
saddle point conditions: 

&=-<£(t-00(*-O> £ 
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(17) 
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Fig. 1: g* /\fkn versus n for different values of <j>. Inset: n c 
versus cf>. For (j> > 0, the continuous curves (corresponding to 
the saddle point that is valid in the region where the unique- 
solution anstaz holds) have been continued as dashed lines 
where multiple flux states satisfy (2) with (3). The highlighted 
y-axis in the inset marks the region in the (<f>, n) phase diagram 
where an expanding phase is possible with g* > occurs. 



7To and ipo denote, respectively, the fraction of "interme- 
diate" metabolites for which the total output equals the 
total input and the fraction of reactions with zero flux: in 
particular, 
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Equations (17) can be solved numerically to obtain g* (and 
all order parameters) as a function of n. Note that the de- 
pendence of g* on the disorder variables is embodied in 
the combination g* /V~k, meaning that the optimal growth 
rate is larger the bigger is the spread in the difference be- 
tween input and output coefficients (in other words, the 
system maximizes growth by taking advantage of imbal- 
ances between inputs and outputs in the stoichiometric 
coefficients). Results for g*/^/nk as a function of n are 
shown in Fig. 1. The critical point n c where g* — (now 
a function of </>) becomes larger as <p increases and ulti- 
mately, when all metabolites form a conserved pool, be- 
comes equal to two. This signals that conserved metabolic 
pools reduce the optimal growth capabilities, e.g. twice as 
many reactions are required to ensure a steady state with 
g = when 0=1 than when = 0. Similarly, the quan- 
tities ip an d tto defined in (20) and (21) are displayed in 
Fig. 2. Notice that as </> increases a larger and larger frac- 
tion of reagents is fully re-cycled into production, while 
the fraction of zero fluxes becomes consistently smaller as 
increases (as more active reactions are required to keep 
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Fig. 2: Fraction of null reactions (ipo, top panel) and fraction of 
intermediate reagents (tiq, bottom panel) versus n for different 
values of 0. Continuous/dashed lines convention as in Fig. 1. 



a steady growth rate). 

These results can be confirmed numerically by calculat- 
ing g* via the MinOver + algorithm introduced in [8] based 
on [11] (see Fig. 3). Its main idea is to rotate the vector 
{si} iteratively in the direction of the worst-satisfied con- 
straint, until all constraints are satisfied at fixed p. Then 
one can increase p by a quantity bp and iterate the pro- 
cedure, until no more solutions are found. Simulations 
suggest that as g — > g* the solution is unique for n < n c 
(where g* < 0), while for n > n c the assumption of a 
unique solution that underlies the x ~~ >• limit ceases to 
be correct and multiple solutions must exist. This can be 
seen directly in numerical experiments by measuring the 
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Fig. 3: Analytical (line) and numerical (markers) results for 
g* versus n. The numerical solution has been calculated by 
MinOver + [8] on systems with NM = 10 4 (variable sizes: 44 < 
N < 216, 46 < M < 220). Corresponding values of are 
displayed in the different panels. Averages over 50 disorder 
samples. The intrinsic resolution of p used for these computer 
experiments is 5p — 10~ 4 . 
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between different solutions s Q and sp obtained by 
MinOver + . If s a = sp, then q a p = 1. In general, 
< q a p < I- (Note that the definition (22) must be cor- 
rected to deal with null fluxes, since the overlap of the null 
solution with itself must obviously equal one.) In particu- 
lar, we are interested in the average of q a p (averaged over 
many solutions), denoted by q. This is reported in Fig. 4. 
It is natural to expect that as p increases, the volume of 
solutions decreases and correspondingly q increases. We 
see that indeed the volume initially shrinks but it stops 
contracting before the optimal growth rate is reached. Ul- 
timately, q < 1 at p* , confirming that many flux states sat- 
isfy (2). By contrast, when no conserved pools are present 
q keeps growing and tends to 1 as p — ► p* , as a unique 
solution exists. 

From a biological viewpoint, one sees that stoichiomet- 
ric constraints on one hand reduce the production capabili- 
ties limiting the system to a stationary state. On the other 
hand, they increase robustness, since many (microscopic) 
flux states are compatible with the optimal (macroscopic) 
growth performance. Clearly, increased flexibility is cru- 
cial for biological stability, since environmental changes or 
localized flux variations are more likely to be sustained by 
the system by re-arranging fluxes while keeping maximal 
production rates. This trade-off between stability and op- 
timal growth has been recently observed in the metabolic 
network of the bacterium E. coli [12]: a finite volume of 
flux vectors is indeed compatible with a maximal growth 
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Fig. 4: Average overlap between 100 different solutions of (2) 
versus p for a network with N = 100 and M = 50 (n — 2) in 
the presence of a single conserved pool formed by 20 metabo- 
lites ((/> — 0.4). For this system, p* = 1. Different lines 
corresponding to different disorder realizations are shown, to 
stress sample-independence of the qualitative picture. Inset: 
same, without the conserved pool (<f> = 0). For this system, 
p* ~ 1.066. 



assumption for E. coli and the predicted flux range is in 
good agreement with the measured experimental fluxes in 
different environments. 

It is easy to extend these results to the case in which 
the reaction network presents dispersion at reaction nodes. 
The simplest possibility is to consider, instead of (3), the 
condition 

M 

Y j z»{a1-b1)=r li (23) 

with rjiS real constants. A positive (respectively, negative) 
r/i says that for reaction i more (respectively, less) reagents 
are produced than they are consumed. While the typical 
properties of large, random instances with random rjiS 
can be characterized analytically along the lines described 
above, it is more useful to concentrate on the case where 
Vi = V which provides for an immediate interpretation 
of the limiting cases ij > and rj < 0. It is easy to 
understand that (23) implies that the average input and 
output stoichiometric coefficients must be linked by 



a — b = 
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This in turn yields 

P*<Pn 
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Fig. 5: Analytical (lines) and numerical (markers) values of 
p* versus n for different r\. The numerical solution is for a 
system of M = 100 metabolites. Analytical lines are given by 
p(n) — p„(l + g* /V nM), with g* given by the solution of (17). 
The horizontal dashed lines are placed at p . 



Hence, the influence of 77 is essentially to allow, in a finite 
system, for an expanding phase if r\ > 0. If r\ < 0, instead, 
the system is necessarily confined to a contracting regime. 
In particular, the introduction of 77 does not affect g*. An- 
alytical and numerical studies fully confirm this prediction 
(see Fig. 5). 

To summarize, the existence of conserved metabolic 
pools (groups of reagents whose aggregate concentration is 
invariant in time) stems directly from stoichiometric bal- 
ance and characterizes all biochemical reaction networks. 
The model we addressed aims at understanding how these 
invariant structures affect the global growth properties 
within a simple exactly solvable setting. We have seen that 
stoichiometry frustrates the system by reducing its opti- 
mal production capability. At the same time, a finite vol- 
ume of flux configurations is compatible with the optimal 
conditions, implying an increase of stability and flexibility. 
Both ingredients are crucial at the biological level since, 
reasonably, metabolic networks should be robust to local 
flux variations. The basic mechanism by which local con- 
servation laws build up to regulate the trade-off between 
growth and stability in the present model is likely to lie at 
the core of recently observed properties of real metabolic 
networks. An important aspect that cannot be analyzed 
theoretically at the level of a fully-connected model is how 
dynamical invariants affect the way in which reaction net- 
works re-arrange fluxes in response to knock-outs or en- 
vironmental changes. It is reasonable to expect that the 
graphical version of the problem, introduced in [8] , will be 



useful to shed some light on this issue. 



We are indebted with S. Franz and T. Jorg for useful 
comments and questions. 
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